#ifndef DATA_STRUCT_CPP
#define DATA_STRUCT_CPP

#include "AlgoEddyDetection.h"
#include <iostream>
#include <Eigen/Dense>
#include <nlohmann/json.hpp>

namespace EddyDetection {

    struct VelocityField {

        std::vector<Eigen::MatrixXd> u_layers; // Each element is a layer of u velocity
        std::vector<Eigen::MatrixXd> v_layers; // Each element is a layer of v velocity
        Eigen::VectorXd lat;
        Eigen::VectorXd lon;
        Eigen::VectorXd depth;
        Eigen::VectorXd time;

        const Eigen::MatrixXd& u() const {
            static Eigen::MatrixXd empty_matrix;
            return u_layers.empty() ? empty_matrix : u_layers[0];
        }
        const Eigen::MatrixXd& v() const {
            static Eigen::MatrixXd empty_matrix;
            return v_layers.empty() ? empty_matrix : v_layers[0];
        }

        VelocityField() {}

        VelocityField(const Eigen::MatrixXd& u_data,
            const Eigen::MatrixXd& v_data,
            const Eigen::VectorXd& lat_data,
            const Eigen::VectorXd& lon_data)
            : lat(lat_data), lon(lon_data) {
            u_layers.push_back(u_data);
            v_layers.push_back(v_data);
        }

        VelocityField(const std::vector<Eigen::MatrixXd>& u_data,
            const std::vector<Eigen::MatrixXd>& v_data,
            const Eigen::VectorXd& lat_data,
            const Eigen::VectorXd& lon_data,
            const Eigen::VectorXd& depth_data = Eigen::VectorXd(),
            const Eigen::VectorXd& time_data = Eigen::VectorXd())
            : u_layers(u_data), v_layers(v_data),
            lat(lat_data), lon(lon_data),
            depth(depth_data), time(time_data) {
        }
    };

    struct TemperatureField {
        std::vector<Eigen::MatrixXd> data_layers;  // Each element is a layer of temperature
        Eigen::VectorXd lat;
        Eigen::VectorXd lon;
        Eigen::VectorXd depth;

        TemperatureField() = default;

        TemperatureField(const Eigen::MatrixXd& temp_data,
            const Eigen::VectorXd& lat_data, const Eigen::VectorXd& lon_data)
            : lat(lat_data), lon(lon_data) {
            data_layers.push_back(temp_data);
        }

        TemperatureField(const std::vector<Eigen::MatrixXd>& temp_data,
            const Eigen::VectorXd& lat_data, const Eigen::VectorXd& lon_data,
            const Eigen::VectorXd& depth_data = Eigen::VectorXd())
            : data_layers(temp_data), lat(lat_data), lon(lon_data), depth(depth_data) {
        }

        const Eigen::MatrixXd& data() const {
            return data_layers.empty() ? Eigen::MatrixXd() : data_layers[0];
        }
    };

    struct EddyContour {
        std::vector<Eigen::Vector2d> boundary_points;  // (lat, lon)
        Eigen::Vector2d center;
        double area_km2;
        double radius;
        double mean_vorticity;
        double mean_omega;
        double temp_anomaly;
        EddyContourDataCyclone type;
        EddyContourDataThermal thermal_type;
        std::string id;
        double depth;
        int layer_index;
    };

    struct Eddy3D {
        std::string id;
        std::string type;
        std::vector<EddyContour> layers;  // A 2D eddy layer that forms a 3D eddy
        int layer_count;
        Eigen::Vector3d center_3d;  // (lat, lon, depth)
        Eigen::Vector2d vertical_extent;  // (min_depth, max_depth)
        double volume_km3;
    };

    namespace Utils {
        // Create test data (dual eddy current)
        VelocityField createDoubleGyreField(int grid_x = 200, int grid_y = 100, double t = 0.0);

        // Data validation
        bool validateVelocityField(const VelocityField& field);
        bool validateTemperatureField(const TemperatureField& field);

        EddyDetectionConfig suggestOptimalParameters(const VelocityField& velocity_field);
        void printDataCharacteristics(const VelocityField& velocity_field);

        // Coordinate System Conversion
        double degToKm(double deg_diff, double mean_lat_rad);
        Eigen::Vector2d pixelToCoordinate(const Eigen::Vector2i& pixel,
            const Eigen::VectorXd& lat,
            const Eigen::VectorXd& lon);
    }

    // 主接口类 - 简化的API
    class EddyDetector {
    public:
        explicit EddyDetector(const EddyDetectionConfig& config = EddyDetectionConfig{});

        std::string detectEddiesLagrangian(const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            double depth_m = 0.0,
            int time_index = 0);

        std::string detectEddiesEulerian(const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            double depth_m = 0.0,
            int time_index = 0);

        std::pair<std::string, std::string> detectEddiesBothMethods(
            const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            double depth_m = 0.0,
            int time_index = 0);

        std::vector<std::vector<EddyContour>> detectLayers2DEddies(const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            const std::string& method = "Eulerian");
        
        std::vector<Eddy3D> detect3DEddies(const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            const std::string& method = "Eulerian");

    private:
        std::vector<EddyContour> detect2DEddiesForLayer(
            const Eigen::MatrixXd& u, const Eigen::MatrixXd& v,
            const Eigen::MatrixXd& temp, const Eigen::VectorXd& lat,
            const Eigen::VectorXd& lon, const std::string& method);

        EddyDetectionConfig config_;
    };

    class DataTransform {
    public:
        void data_to_field(const VelocityData& velocity_data, const TemperatureData& temp_data,
            VelocityField& velocity_field, TemperatureField& temp_field);
        void eddy2D_to_2DData(std::vector<std::vector<EddyContourData>>& eddiesdata_3d, const std::vector<std::vector<EddyContour>>& eddies_3d);
        void eddies_to_eddiesdata(std::vector<Eddy3DData>& eddiesdata_3d, const std::vector<Eddy3D>& eddies_3d);
    };

    class LagrangianCore {
    public:
        static Eigen::MatrixXd computeFTLE(const VelocityField& velocity_field,
            const EddyDetectionConfig& config,
            int time_index = 0);

        static std::vector<EddyContour> detectEddiesLagrangian(
            const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            const EddyDetectionConfig& config,
            int time_index = 0);

    private:

        static Eigen::MatrixXd computeFlowMap(const VelocityField& velocity_field,
            double t_start, double t_end,
            const EddyDetectionConfig& config);

        static std::vector<std::vector<Eigen::Vector2d>> findFTLEContours(
            const Eigen::MatrixXd& ftle_field,
            const Eigen::VectorXd& lat, const Eigen::VectorXd& lon,
            double threshold_percentile);
    };

    class EulerianCore {
    public:

        static Eigen::MatrixXd computeVorticity(const VelocityField& velocity_field,
            int time_index = 0);


        static std::vector<EddyContour> detectEddiesEulerian(
            const VelocityField& velocity_field,
            const TemperatureField& temp_field,
            const EddyDetectionConfig& config,
            int time_index = 0);

    private:
        static std::vector<std::vector<Eigen::Vector2d>> findVorticityContours(
            const Eigen::MatrixXd& vorticity_field,
            const Eigen::VectorXd& lat, const Eigen::VectorXd& lon,
            double threshold_percentile, double min_area_km2);

        static std::vector<std::vector<Eigen::Vector2d>> findVorticityContours1(
            const Eigen::MatrixXd& vorticity_field,
            const Eigen::VectorXd& lat, const Eigen::VectorXd& lon,
            double threshold_percentile, double min_area_km2);
    };

    class EddyAnalyzer {
    public:

        static EddyContour analyzeEddy(const std::vector<Eigen::Vector2d>& contour,
            const Eigen::MatrixXd& vorticity_field,
            const Eigen::MatrixXd& temp_field,
            const Eigen::VectorXd& lat,
            const Eigen::VectorXd& lon,
            double min_area_km2);


        static std::vector<EddyContour> filterAndSortEddies(
            std::vector<EddyContour>& eddies, double min_area_km2, const VelocityField& velocity_field);

    private:
        static double calculateArea(const std::vector<Eigen::Vector2d>& contour);
        static double calculateRadius(const std::vector<Eigen::Vector2d>& contour);
        static Eigen::Vector2d calculateCenter(const std::vector<Eigen::Vector2d>& contour);
        static double calculateMeanProperty(const std::vector<Eigen::Vector2d>& contour,
            const Eigen::MatrixXd& field,
            const Eigen::VectorXd& lat,
            const Eigen::VectorXd& lon);
    };

    class ThreeDLinker {
    public:
        static std::vector<Eddy3D> link3DEddies(
            const std::vector<std::vector<EddyContour>>& all_layers_eddies,
            const EddyDetectionConfig& config);
    private:
        static Eddy3D create3DEddy(const std::vector<const EddyContour*>& chain);
    };

    class GeoJSONExporter {
    public:

        static std::string exportToGeoJSON(const std::vector<EddyContour>& eddies,
            const std::string& method_name,
            double depth_m);


        static nlohmann::json createGeoJSONObject(const std::vector<EddyContour>& eddies,
            const std::string& method_name,
            double depth_m);


        static std::string export3DEddiesToGeoJSON(
            const std::vector<Eddy3D>& eddies_3d,
            double depth_m);

    private:
        static nlohmann::json createFeature(const EddyContour& eddy,
            const std::string& method_name,
            double depth_m);
        static nlohmann::json create3DEddyFeature(const Eddy3D& eddy_3d);
    };

    static void printVelocityFieldInfo(const VelocityField& field) {
        std::cout << "=== 速度场信息 ===" << std::endl;
        std::cout << "纬度范围: " << field.lat.minCoeff() << " 到 " << field.lat.maxCoeff()
            << " (共 " << field.lat.size() << " 个点)" << std::endl;
        std::cout << "经度范围: " << field.lon.minCoeff() << " 到 " << field.lon.maxCoeff()
            << " (共 " << field.lon.size() << " 个点)" << std::endl;

        if (!field.u_layers.empty()) {
            const auto& u = field.u_layers[0];
            const auto& v = field.v_layers[0];
            std::cout << "U速度范围: " << u.minCoeff() << " 到 " << u.maxCoeff() << " m/s" << std::endl;
            std::cout << "V速度范围: " << v.minCoeff() << " 到 " << v.maxCoeff() << " m/s" << std::endl;
        }
        else {
            std::cout << "The velocity field data is empty" << std::endl;
        }

        std::cout << "时间步数: " << field.time.size() << std::endl;
        std::cout << "深度层数: " << field.u_layers.size() << std::endl;
    }

    static void printTemperatureFieldInfo(const TemperatureField& field) {
        std::cout << "=== 温度场信息 ===" << std::endl;
        std::cout << "纬度范围: " << field.lat.minCoeff() << " 到 " << field.lat.maxCoeff()
            << " (共 " << field.lat.size() << " 个点)" << std::endl;
        std::cout << "经度范围: " << field.lon.minCoeff() << " 到 " << field.lon.maxCoeff()
            << " (共 " << field.lon.size() << " 个点)" << std::endl;

        if (!field.data_layers.empty()) {
            const auto& data = field.data_layers[0];
            std::cout << "温度范围: " << data.minCoeff() << " 到 " << data.maxCoeff() << " °C" << std::endl;
        }
        else {
            std::cout << "The temperature field data is empty" << std::endl;
        }

        std::cout << "深度层数: " << field.data_layers.size() << std::endl;
    }

    static void printEddyInfo(const EddyContour& eddy) {
        std::cout << "=== Eddy information ===" << std::endl;
        std::cout << "ID: " << eddy.id << std::endl;
        std::cout << "Type: " << eddy.type << std::endl;
        std::cout << "Thermodynamics type: " << eddy.thermal_type << std::endl;
        std::cout << "Central coordinates: (" << eddy.center[0] << ", " << eddy.center[1] << ")" << std::endl;
        std::cout << "Area: " << eddy.area_km2 << " km²" << std::endl;
        std::cout << "Average vorticity: " << eddy.mean_vorticity << "/s" << std::endl;
        std::cout << "Average omega: " << eddy.mean_omega << "/s" << std::endl;
        std::cout << "Abnormal temperature: " << eddy.temp_anomaly << " °C" << std::endl;
        std::cout << "Boundary point count: " << eddy.boundary_points.size() << std::endl;
    }

    static void printConfigInfo(const EddyDetectionConfig& config) {
        std::cout << "=== Detecting configuration information ===" << std::endl;
        std::cout << "FTLE Threshold percentile: " << config.ftle_threshold_percentile << "%" << std::endl;
        std::cout << "Vorticity Threshold percentile: " << config.vorticity_threshold_percentile << "%" << std::endl;
        std::cout << "Minimum area threshold: " << config.min_area_km2 << " km²" << std::endl;
        std::cout << "Integration time: " << config.integration_time_hours << " hour" << std::endl;
        std::cout << "Use OpenMP: " << (config.use_openmp ? "Yes" : "No") << std::endl;
        std::cout << "Number of threads: " << config.num_threads << std::endl;
    }

} // namespace EddyDetection 

#endif // !DATA_STRUCT_CPP
